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Abstract: In gold-gold collisions of the Relativistic Heavy Ion Collider a perfect fluid, the strongly interacting quark 

gluon plasma (sQGP) is created. The time evolution of this fluid can be described by hydrodynamical 
models. After an expansion, hadrons are created during the freeze-out. Their distribution reveals informa- 
tion about the final state. To investigate the time evolution one needs to analyze penetrating probes: e.g. 
direct photon observables. In this paper we analyze a 1-1-3 dimensional solution of relativistic hydrody- 
namics. We calculate momentum distribution, azimuthal asymmetry and momentum correlations of direct 
photons. Based on earlier fits to hadronic spectra, we compare photon calculations to measurements to 
determine the equation of state and the initial temperature of sQGP. We find that the initial temperature 
in the center of the fireball is 507±12 MeV, while for the sound speed we get Cs=0.36±0.02. We also 
estimate a systematic error of these results. We find that the measured azimuthal asymmetry is also not 
incompatible with this model, and predict a photon source that is significantly larger in the out direction 
than in the side direction. 
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1. Introduction 

In the last several years it has been revealed that the matter produced in the collisions of the Relativistic Heavy 
Ion Collider (RHIC) is a nearly perfect fluid [1], i.e. it can be described with perfect fluid hydrodynamics. There 
was a long search for exact hydrodynamic models (solutions of the partial differential equations of hydrodynamics) 
and several models proved to be applicable [2-10]. Many solve the equations of hydrodynamics numerically, which 
has the advantage of having arbitrary initial conditions. However, it is possible to find analytic solutions as well, 
see the above cited references for example. These solutions will have a couple of parameters, the optimal value of 
which can then be determined by fitting theoretical results to data. Among these models there are non-relativistic 
l-f3 dimensional ones, as well as 1-fl dimensional relativistic models - but not many l-f3 dimensional and 
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relativistic models were tested yet. In this paper we extract photon observables from the relativistic, ellipsoidally 
symmetric model of Ref. [11]. Hadronic observables were calculated in Ref. [12]. Here we calculate transverse 
momentum distribution and azimuthal asymmetry (elliptic flow) of direct photons. We also calculate Bose- 
Einstein correlations and their widths, called HBT radii. 

For the direct photon calculations, the framework used in this paper has to be made very clear. Our key assumption 
is, that even though direct photons may not be thermalized in the strongly interacting plasma (as their mean free 
path may be on the order of the size of the fireball), but the radiation itself is thermal. Thus the phase-space 
distribution of the photons is characterized by the temperature of the medium (at a given space-time cell), while 
the expansion of the fireball also effects the observed spectrum. This is a macroscopic model, and in the following 
we will calculate photon observables from it and compare to RHIC data. The most important assumption is, that 
the spectrum of direct photons is thermal because macroscopically, the photon radiation is thermal. Our tool will 
be an exact solution of relativistic hydrodynamics, applied to the expanding medium, the strongly interacting 
plasma created in Au+Au collisions at RHIC. 

2. Perfect fluid hydrodynamics 

Perfect fluid hydrodynamics is based on local conservation of a conserved charge (n) and energy- momentum (T*^"): 

a^(nu'^) = 0, (1) 
^mT""^ = 0, (2) 

where is the flow field in the fluid. These quantities here refer to the strongly interacting plasma created at 
RHIC, i.e. the equations dealing with its temperature, flow etc. The fluid is perfect if the energy-momentum 
tensor is diagonal in the local rest frame, i.e. viscosity and heat conduction are negligible. Thus T'"' is chosen as 

T"" = (e + p)^t''u'^-pff''", (3) 

where e is energy density, p is pressure and g''" is the metric tensor, diag(-l, 1,1,1). Note also that in this 
paper, x*^ = {t,rx-,ry,r^) is a given point in space-time, 9^ = d/dx'^ is the derivative versus space time, while 
p'^ — {E,px,Py,p~) is the four-momentum. 

The conservation equations are closed by the equation of state (EoS), which gives the relationship between energy 
density e and pressure p. Typically e = up is chosen, where the proportionality "constant" k may depend on 
temperature T. Note also that k equals l/Cj, with Cs being the speed of sound. Temperature in turn is connected 
to density n and pressure p via p — nT. In some solutions (as also on the analyzed one) a bag constant B can be 
introduced (this is however not favored by the data, because no first order phase transition is seen in high energy 
heavy ion collisions, these are in the cross-over regime of the QCD phase-map [13]). The exact, analytic result 
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for hydrodynamic models is, that the hadronic observables do not depend on the initial state or the dynamical 
equations separately, just through the final state [14]. But if we fix the final state from hadronic data, we can 
determine initial state parameters from direct photon spectra. 

Even though many solve the above equations numerically, there are only a few exact solutions for these equa- 
tions. One (and historically the first) is the implicit solution discovered more than 50 years ago by Landau 
and Khalatnikov [2-4]. This is a 1+1 dimensional solution, and has realistic properties: it describes a 1+1 di- 
mensional expansion, does not lack acceleration and predicts an approximately Gaussian rapidity distribution. 
Another renowned solution of relativistic hydrodynamics was found by Hwa and Bjorken [5-7]: it is simple, 1+1 
dimensional, explicit and exact, but accelerationless. It is boost-invariant in its original form. Boost invariance 
is however incompatible with data from RHIC, so the solution fails to describe the data, it still can be used to 
estimate the energy density in high energy heavy ion collisions. Important are solutions [9, 10] which are explicit 
and describe a relativistic acceleration, i.e. combine the properties of the Landau-Khalatnikow and the Hwa- 
Bjorken models. With these models one can have an advanced estimate on the energy density, but investigation 
of transverse dynamics is not possible. 

There was only one 1+3 dimensional relativistic solution investigated: the solution in Ref. [11]. Hadronic ob- 
servables from this solution were computed and compared to data in Ref. [12]. Present paper calculates thermal 
photon observables from this realistic 1+3 dimensional model and compares them to data for the first time. Our 
method is different from numerical calculations: here one can determine the best values of the parameters of the 
solution by fitting the analytic model results to data. 

3. The analyzed solution 

The analyzed solution [11] assumes self-similarity and ellipsoidal symmetry, as described in Ref. [12]. The ellip- 
soidal symmetry means that at a given proper time the thermodynamical quantities are constant on the surface 
of expanding ellipsoids. The ellipsoids are given by constant values of the scale variable s, defined as 

^2 2 2 

+ + (4) 



x(t)^ y(t)2 z{ty' 

where X{t), Y{t), and Z{t) are time dependent scale parameters (axes of the s = 1 ellipsoid), only depending on 
the time. Spatial coordinates are rx, ry, and r^. The velocity-field (of the fireball created in RHIC collisions) is 
described by a Hubble-type expansion: 

u'^M-yfl mr Mr] (5) 

U [Xj-J ^1, ^^^^ Tx, y^^^ Ty, ^^^^ T^j , (5) 

where x means the four-vector {t,rx,ry,rz), and X{t) = dX{t)/dt, similarly for Y and Z. The X{t) — Xq, 
Y{t) = Xo, Z{t) = Xo (i.e. all are constant) criteria must be fulfilled, ie. the solution is accelerationless. This is 
a drawback of this model. 
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The temperature distribution T{x) is 

where r is the proper time, s is the above scaUng variable, !/(s) is an arbitrary function, while To = T\s=o,t^to, 
and To is an arbitrary proper time, but practically we choose it to be the time of the freeze-out, thus To is the 
central freeze-out temperature. The function u{s) is chosen as u{s) = exp(— 6s/2) where b is proportional to the 
temperature gradient. If the fireball is the hottest in the center, then 6 < 0. If there is a conserved charge in 
the system e.g. the baryon number density, then charge number density n{x) can be utilized in the solution. As 
described in Refs. [11, 12], such a number density can be introduced as 



n(x) = no (J^y u{s). (7) 



For the momentum distribution of direct photons, this will not be needed, as the only the temperature of the 
medium (the strongly interacting plasma) governs the creation of photons, not the density (which however plays 
an important role also in the case of hadron creation). 

4. Hadronic and photonic source functions 

The picture used in hydro models is that the pre freeze-out (FO) medium is described by hydrodynamics, and the 
post freeze-out medium is that of observed hadrons. The hadronic observables can be extracted from the solution 
via the phase-space distribution at the FO. This will correspond to the hadronic final state or source distribution 
S{x,p). See details about this topic in Ref. [12]. It is important to see that the same final state can be achieved 
with different equations of state or initial conditions [14]. However, as discussed below, the source function of 
photons is sensitive to the whole time evolution, thus both to initial conditions and equation of state as well. 
In our model the hadronic source distribution of bosons takes the following form [12]: 

S{x,p)d X = M , , , ; , ' -, (8) 

^ ' n{x)exp{p^ui'{x)/T{x)) -1 

where JV — g/{2Tv)'^ (with g being the degeneracy factor), H{t) is the proper-time probability distribution of the 
FO. It is assumed to be a 5 function or a narrow Gaussian centered at the freeze-out proper-time tq. Furthermore, 
IJ,{x)/T{x) — lnn(a;) is the fugacity factor and d^E^(a-)p'' is the Cooper-Frye factor [15] describing the fiux of the 
particles, and d'^'S^{x) is the vector-measure of the FO hyper- surface. Here the source distribution is normalized 
such as J S{x,p)d'^xd'^p/ E = A'^, i.e. one gets the total number of particles N (using c=l, h=l units). 
For the source function of photon creation we have 

S{x, p)d^x = M PMrf^S^Wrf^ ^ ^ P^^U^ ^4^ 

^ '-^^ e^p{p^ut^{x)/T{x))-l ex.p{p^,ut^{x)/T{x))~l 
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where p^d'^E'' is again the Cooper-Frye factor of the emmission hypersurfaces. Similarly to Ref. [12] we assume 
that the hyper-surfaces are parallel to u'^ , thus d''E^(a;) = u^dPx. This yields then P/ju'' which is the energy 
of the photon in the co-moving system. Here the simple assumption of a thermal distribution is used, i.e. that 
the radiation of thermal photons follows the usual oc E/ [exp(_E/r) — 1) type of behavior. Photon production 
time has to be also taken into account, and we assume that photons are created from an initial time ti until a 
point sufficiently near the quark-hadron transition. Note however, that in this case the photon production rate 
(integrated over momentum) is not proportional to T^, as expected from some microscopic models, where 



with ua and ub being the densities of the two input particles, while (ja+b^x is the cross-section of the given 



give a production rate proportional to T , or even larger if the velocity for example rises with temperature. This 
wil bel incorporated in our model as a multiplicative factor in the photon source function, but we will use this 
only to estimate a systematic uncertainty of the parameters of our model when comparing to real data. Note that 
a similar argument is used in UrQMD, where all possible microscopic processes are taken into account, but the 
measured pion, kaon and proton spectra do not have such a good agreement with this model as the agreement 
with hydrodynamic models. 

Experimental observables can then be calculated from the source function, most importantly the invariant mo- 
mentum distribution Ni (p) as 



where again the integration in time is on a definite interval, from the initial time until tq. To perform this inte- 
gration we use a second order saddlepoint approximation. In this approximation the point of maximal emissivity 
is 



rate(A + B ^ X) = uaub {ga+b^xv) , 



(10) 



process and v is the incoming velocity. Here, however, the densities are proportional to T"^ normally, so this would 




(11) 




(12) 



ro.v = 




(13) 




(14) 



while the widths of the particle emitting source are 




(15) 



(16) 



(17) 
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where we introduced the auxihary quantities 

Py = (19) 

'o 

P. = (20) 

where again k = cj^ is describing the EoS. Note that the source width depend clearly on time, as the system is 
expanding and so does the photon source. 

5. Thermal photon observables 

The invariant one-particle momentum distribution is defined in eq. (11), it depends on the three-momentum 
V = {PTjPy,Pz)- We will introduce the {pt, ^p,Pz) cylindrical coordinates {z being the beam direction) and use the 
longitudinal rapidity y (for which E dy — dp^ is true). The PHENIX detector is aligned such that \y\ < 0.35 for 
photons. Thus, as usual, we restrict our calculations to y = from now on (note that in this case E — pt is true 
for photons). Our key quantity is then 

iV.(p)U-A^.(P.^) = ^ (21) 

We transform this into two one-dimensional observables, the invariant transverse momentum distribution and the 
elliptic flow, similarly to Ref. [12]: 

Ni{pt, ip) = Ni{pt) ^1 + 2 v„ cos(nv5) j (22) 
where V2, the second Fourier component, is the elliptic flow. Ni{pt) and V2(pt) can be calculated from Ni{pt,'p) 



Niipt) = / Ni{p)dip (23) 
Jo 

^ f^''Ni(p)cos(2ip)dip 

V2{pt) = ^ , ; (24) 

Ni(pt) 

In order to calculate these quantities, we integrate on the space-time coordinates. We also have to integrate on 
the azimuthal angle (p. During that, similarly to Ref. [12] we have to use the modifled Bessel functions /„ The 
result on the invariant transverse momentum distribution depends on the initial and flnal times. We introduce 
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the variable ^ = (with to as the freeze-out time), then in terms of ^ the time integration goes from i to 1, if 



Ni{Pt) = ^(27r) VP-PyP^^oTo 



To 



3yl"+¥-i 



(Cao„ + Dai,0 r ( n + — - g ' ^fr?" 



(25) 



where we introduce 



A = 1 



D 



px + I 



B = 



Px - Py 



Q + Py _|_ P^ + pI 



2 _ 2 
Py _|_ Pa^ Py 



(26) 
(27) 
(28) 
(29) 



and ao, ai are the Taylor-coefficients of the first two modified Bessel functions: 



""=(^'0'i'°'^'0'23^'°' 
(°'^'°']^'°'3^'°'T8i32' 



(30) 
(31) 



The modified Bessel functions can be expressed then as 



Io{x) = y^^aonx" 

n = 
oo 



(32) 
(33) 



As the coefficients are strongly decreasing, in real calculations we can restrict ourselves to use only the first two 
of them, i.e. we can make the approximation of Io{x) = x and /i(a:) = x^/2. 
For the elliptic fiow we get: 



1 (2vr) 2 ^pxPyPz ^ ( Pt\^ 



(34) 



C-3, , C-1 

— - — (aon + a2n) H — ain 

4 2 



rl^-L^'^ 3 .Pt . Pr+Py + P£ ,4k 5 .Pt.3 
i nH ,A— f " H ain^l nH 



Here /2 is the modified Bessel function, which has the Taylor-coefficients of a2 — (O, 0, ^,0, ^,0, ^^2^^^ ■■•)• 
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Parameter 

hline Central FO temperature To 



FO proper-time 
Eccentricity 
Transverse expansion 
Longitudinal expansion 
Equation of State 
Initial time 



TO 
€ 

zl/b 



Value 
204 MeV 
7.7 fm/c 
0.34 
-0.34 
-1.6 
7.9 ±0.7 



Type 
fixed 
fixed 
fixed 
fixed 
fixed 
free 



ti — 0.7 fm/c interval 



Fit property 



Value 



Number of data points 
Fitted parameters 
Degrees of freedom 
Chisquare 
Confidence level 



5 
2 

NDF 5 - 2 : 
7.0 
7.2% 



Table 1. Parameters of the solution, describing the expanding sQGP. The first five are taken from the hadronic fits of Ref. [12] 
(see more details about the parameters in the reference). The EoS parameter k is fitted, while for the initial time we 
determine an interval of acceptability (with 95% confidence). 

6. Comparison to the measured direct photon spectrum 

Hadronic data were already described with this model in Ref. [12]. Freeze-out properties were thus determined 
from hadronic fits. These properties include the expansion rates, the freeze-out proper-time and freeze-out 
temperature (in the center of the fireball), see Table 1 for the values of these parameters. We use the parameters 
of the hadronic fit and leave only the remaining as free parameters. The free parameters will be At (the equation 
of state parameter) and ti, the initial time of the evolution. 

Similarly to Ref. [12], we use transverse expansion {ut) and eccentricity (e) instead of x and y direction expansion 
rates: 

1 _ 1 / 1 1 

We use direct photon data from the PHENIX Collaboration [16], measured in ^/snn =200 GeV Au+Au collisions. 
For the fits we use the standard Minuit package [17]. Note that we also utilize a normalizing factor to describe 
the data. Fit parameters are summarized in Table 1. The fit itself is shown on Fig. 1. 
The equation of state result is «: = 7.9 ± 0.7, or alternatively, using k = l/c^: 

Cs = 0.36 ± 0.02 (36) 

which is in nice agreement with both lattice QCD calculations [18] and experimental results from hadronic 
data [19, 20]. This represents an average EoS as it may vary with temperature. There may be solutions with a 
K(r) function, but for the sake of simplicity we assumed here an average, fixed Note however that the spectrum 
is not very sensitive to the initial time as in early times the thermal photon emission is not in the region of the 



X2 + r2 ■ 



(35) 
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2 2.5 

Pt (GeV/c) 



Fi gure 1. Invariant transverse momentum of direct photons from our hydro model, normalized to the data. The model validity 
goes until roughly 3 GcV in transverse momentum, the last point on the plot is included in the fit to have enough 
constraints. 
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Figure 2. Invariant transverse momentum of direct pliotons, calculated with difl^erent initial times. This figure illustrates the 
inscnsitivity on the initial time. 



experimental data. See Fig. 2 as an illustration: we plotted the invariant transverse momentum distribution 
Ni{pt) for different values of initial time. Clearly the very first times do not have a large contribution to N\{pt) 
in the desired region. Our model does not contain acceleration, it is a Hubble-fiow type of model, but the initial 
acceleration does not play a large role in the thermal photon spectrum because of this insensitivity on the initial 
time. As a result, the fits will neither be sensitive to the exact value of the initial time - or otherwise, the initial 
time ti will have a very large error. As a workaround, we determine an "interval of acceptability" for ti. The 
maximum value for ti within 95% probability is 0.7 fm/c. This can then be used to determine a lower bound for 
the initial temperature, using the eq. 6. Thus the initial temperature of the fireball (in its center) is: 



Ti > 507 ± 12MeV 



(37) 
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at 0.7 fm/c. The uncertainty comes from the uncertainty of k. This is in accordance with other hydro models as 
those values are in the 300 — 600 MeV interval [16] . 

As a systematic test, we investigate the changes in the result when using a factor of (T/To)^, with iV = 0, 1, 2, 3. 
Such a factor would arise if the photon creation can be described by a microscopic process mentioned in Section 4. 
Such a prefactor causes only minor change in the resulting spectrum, as it is dominated by the exponential factors 
in it. However, the equation of state parameter k changes from 7.9 to 6.5 as we increase the exponent in the 
prefactor. It is important also, that the best fit (smallest x^) is achieved without the prefactor. This test can still 
be used as a systematic error of our estimate, thus the final result for Cs is then 

Cs = 0.36 ± 0.02,t„t ± OMsyst (38) 

while for the temperature before 0.7 fm/c we get 

Ti > 507 ± Ustat ± 90s,;.tMeV (39) 



7. Direct photon elliptic flow and correlation radii 

A measurement of direct photon elliptic flow was also performed recently at PHENIX [21]. Using the previously 
determined fit parameters we can calculate the elliptic flow of direct photons in Au+Au collisions at RHIC. We 
compare the elliptic flow V2 from eq. 34 to data of Ref. [21], using the parameters of Table 1, with the exception 
that the eccentricity parameter e has to be modified, as the hadronic elliptic flow data of Ref. [12] utilized also a 
different e. Due to the low number of points in the desired range, a fit can not be performed here, but we use the 
average value e in case of the two fits of Ref. [12]. The resulting curve, where the value e = 0.59 is used, is shown 
on Fig. 3. Note that already for 1 GeV/c the elliptic fiow is very small. It is quite unusual, but V2 is negative 
for even smaller transverse momenta. This is due to the fact that in the photon source, the Pi_iu'^ prefactor (the 
comoving energy) also oscillates in the azimuthal angle tp, as shown in eq. 34. After spatial integration the (p 
dependence will look like: 

A^i (0) cx (1 - Q cos{2(l,))e^ (40) 



where 



(p. - if + [py - if 



{po. - If + {Py -If +2 + 2{pi + pI + pi)/C 



P^iPo: 



Py)\ 



(41) 
(42) 



3/k 



with C = ^ \^) ' "^^^^ function changes sign in the oscillation at a given point, depending on In fact the 
number of oscillations within the [0, 2n] interval will be four (in contrast to two oscillations of a regular cos(2(^) 
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Figure 3. Thermal photon elliptic flow with parameters from Table 1, compared to PHENIX data of Rcf. [21], The eccentricity 
parameter e is talcen from hadronic fits of Ref. [12] as detailed in the text. 



function) if— l<l//3— l/a<l. The two boundaries of this interval can be calculated by solving P = a/{l + a) 
and /3 — a/{l ~ a) for pt (which appears in i^). Exact shape of the elliptic flow depends strongly on input 
parameters, and also it depends stronger on initial time than the invariant momentum distribution. Thus an 
accelerating model might result in a different elliptic flow. 

We also calculated Bose-Einstein correlation radii from the above model. As usual, the two-particle correlation 
function for identical particles can be calculated from the single particle source function S as 



Ciiq) = 1 + A 



~m\ (43) 

5(0) 



where q is the momentum difference of the two particles and S{q) is the Fourier-transformed of the source S{x). 
This transformation can however not be done analytically, so we just calculate the correlation function for the 
parameter set used to describe the direct photon spectrum. We analyze the correlation functions to determine 
the HBT type of correlation radii. We flnd the numerically calculated correlation function to have a peak of 
shape C2(q) = exp where is the matrix of correlation radii, and for a = 2, this represent a Gaussian 

shape. From this shape, we determine the HBT radii -Rout and -Rsido for different average pt values. 
In case of hadronic HBT, correlation radii in the side and out directions are almost equal, as for the hadronic 
transition is of cross-over type (i.e. the transition time is short), see details for example in Ref. [12]. However, in 
case of photons, the creation spans the whole evolution of the flreball, thus i?out will be significantly larger than 
inside- Indeed this is observed in our model, as shown on Fig. 4. 



8. Summary 

The medium created in relativistic Au+Au collisions at RHIC turned out to be a strongly interacting, perfect 
fluid. Observed hadrons are created at the freeze-out of this fluid, while thermal photons are constantly emitted. 
Hadrons thus reveal information about the final state, whereas thermally radiated photons carry information 
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Figure 4. Thermal photon HBT radii (correlation widths) with example parameters from Table 1. Due to the length of photon 
production, i?out is significantly larger than R^idc- Shape of the correlation functions is nearly Gaussian, optimal values 
of a are in the range of 1.7 to 1.9. 

about the whole time evolution. We used direct photon data of PHENIX to determine the equation of state and 
the initial temperature of quark gluon plasma. We find that thermal radiation is consistent with the data, and 
our result on equation of state (average speed of sound) is — 0.36 ± 0.02stt,t ± 0.04syst, and we set a lower 
bound on the initial temperature of the sQGP to 507± 12stat iQOsyst MeV at 0.7 fm/c. To our knowledge, this is 
the first time when these values are extracted from photon and hadron data simultaneously. We also investigated 
the thermal photon elliptic flow from this model, which is not incompatible with measurements. We calculated 
photon HBT radii from the model, and discovered a significantly larger -Rout than i?sido. 
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